#include "header.h"

static fftwq_plan plan_dst, plan_dct;
static __float128 *k, overN, q;
static __float128 coeff, *ra[2], *ro[2];
static __float128 ip;
static int ne, no;

int main(int argc, char** argv)
{
  work	  wrk;
  params  ctrl;
  aux	  data;
  int	  Emm, K, Jay, l;
  int 	  nref, d = 0, n = 0;
  __float128 *alph, *beta, *u, *X, *Y;

  char str[480], line[480], value[480], param[480], value2[480], *dummy;
  sprintf(str, "%s.cfg", argv[1]);

  init_restore(&argc, argv, &ctrl, &wrk);
  FILE* fh = fopen(str,"r"); ctrl.d = 0;
  if (fh == NULL) err_msg("Cannot open file: input.cfg");
  ctrl.n = 0; ctrl.d = 0;
  while (fgets(line, 480, fh)!=NULL) {
    sscanf(line, "%s\t%s", param, value);
    if (strcmp(param,"#n_poles=") == 0) ctrl.d = atoi(value);
    if (strcmp(param,"#nsurfac=") == 0) ctrl.n = atoi(value);
    if (strcmp(param,"#polname=") == 0) sprintf(ctrl.resname,"%s", value);
  }
  fclose(fh);
  d = ctrl.d; n = ctrl.n;

  alph = fftwq_malloc(sizeof(__float128)*d);
  beta = fftwq_malloc(sizeof(__float128)*d);

  u = fftwq_malloc(sizeof(__float128)*n);
  Y = fftwq_malloc(sizeof(__float128)*n);
  X = fftwq_malloc(sizeof(__float128)*n);

  FILE* fhdat = fopen(ctrl.resname,"r");   
  if (fhdat == NULL) err_msg("Cannot open file: poles data file");

  int jj = 0, ind;
  __float128 y0 = 0.Q;
  while (fgets(line, 480, fhdat)!=NULL) {
    if (jj == 0) {
      dummy = fgets(line, 160, fhdat);
      sscanf(line, "# N = %s\tL2_residual = %s\ty0 = %s", param, value, value2); 
      printf("%s\n%s\t%s\t%s\n", line, param, value, value2);
      y0 = strtoflt128(value2, NULL);
      dummy = fgets(line, 160, fhdat);
    }
    if (jj > 0) {
      sscanf(line, "%s\t%s\t%s\t%s\n", param, value, value2, value2);
      ind = jj-1;
      alph[ind] = strtoflt128(value, NULL);
      beta[ind] = strtoflt128(param, NULL);
      printf(" %3d   %.32Qe   %.32Qe\n", ind+1, alph[ind], beta[ind]);  
    }
    jj += 1;
    if (jj == d+1) break;
  }
  fclose(fhdat);

  FILE *fhprf = fopen("profile.dat","w");
  fprintf(fhprf, "# 1. u 2. x-tilde 3. y\n# N = %d\n\n", n);
  for (int jj = 0; jj < n; jj++) {
    u[jj] = 2.Q*pi*jj/n;
    X[jj] = 0.Q;
    Y[jj] = y0;
    for (int dd = 0; dd < d; dd++) {
      X[jj] += crealq(alph[dd]/(tanq(0.5Q*(u[jj] - pi)) - 1.IQ*beta[dd] )); 
      Y[jj] += cimagq(alph[dd]/(tanq(0.5Q*(u[jj] - pi)) - 1.IQ*beta[dd] ));
    }
    fprintf(fhprf, "%.32Qe\t%.32Qe\t%.32Qe\n", u[jj], X[jj], Y[jj]);
  }
  fclose(fhprf);
  
//------------------------------------------------
  int ne, no;
  ne = n/2 + 1; no = n/2 - 1;

  ra[0] = fftwq_malloc(sizeof(__float128)*ne);
  ra[1] = fftwq_malloc(sizeof(__float128)*ne);

  ro[0] = fftwq_malloc(sizeof(__float128)*no);
  ro[1] = fftwq_malloc(sizeof(__float128)*no);

  wrk.XX  = fftwq_malloc(sizeof(__float128)*no);
  wrk.XXk = fftwq_malloc(sizeof(__float128)*no);
  wrk.YY  = fftwq_malloc(sizeof(__float128)*ne);
  wrk.YYk = fftwq_malloc(sizeof(__float128)*ne);

  plan_dst = fftwq_plan_r2r_1d(no, ro[1], ro[0], FFTW_RODFT00, FFTW_ESTIMATE);
  plan_dct = fftwq_plan_r2r_1d(ne, ra[1], ra[0], FFTW_REDFT00, FFTW_ESTIMATE);

  for (int jj = 0; jj < ne; jj++) wrk.YY[jj] = 1.Q*(Y[jj])/n;  
  ifft_even(wrk.YY, wrk.YYk);

  //for (int j = 0; j < no; j++) wrk->XXk[j] = wrk->YYk[j+1];  
  //ifft_odd(wrk->XXk, wrk->XX);  // seg_fault here
  
//--------------------------------------

  init_prep(&ctrl);
  writer_out(wrk.YYk, "restart.rcn", 1., y0);







}
